Cell type classification in tumor microenvironment using HiTME

Background

In this vignette we will use HiTME to classify the cell types found in whole-tumor samples from Breast cancer. This single-cell data set was processed from original work from Bassez et al. (2021) Nature Medicine. The original data set contains single-cell data from hormone receptor-positive or triple-negative breast carcinoma biopsies from 42 patients. One cohort (n=31) of patients with non-metastatic, treatment-naive primary invasive carcinoma of the breast was treated with one dose of pembrolizumab (Keytruda or anti-PD1) approximately 9 ± 2 days before surgery. A second cohort of patients (n=11) received neoadjuvant chemotherapy for 20–24 weeks, which was followed by anti-PD1 treatment before surgery. In both cohorts, a tumor biopsy was harvested immediately before anti-PD1 treatment (‘pre’), while another biopsy was collected during subsequent surgery (‘on’).

Full Data Original Source

Demo dataset link

Goals

Use HiTME tools to predict cell types in 3 different whole tumor samples of breast cancer.
Evalute their fitness and compare to authors prediction

R Environment

# install HiTME if needed

# remotes::install_github('carmonalab/HiTME')

# load libraries
library(HiTME)
library(Seurat)

# increase timeout max when downloading dataset and references
options(timeout = max(2000, getOption("timeout")))

Processed dataset

A Seurat object with 3 patient samples from pre-treatment (‘pre’) condition were subset from the original data set Bassez et al. 2021: “BIOKEY_13_Pre”, “BIOKEY_14_Pre”, and “BIOKEY_5_Pre”, without downsampling or further per sample filtering. No cells were discarded based on quality control metrics and data was normalized using Seurat NormalizeData function. All available metadata was included.

File can be downloaded from (https://figshare.com/ndownloader/files/43848153)

Demo data downloading

ddir <- "input"
dest <- file.path(ddir, "bassez_3patients.rds")

if (!file.exists(dest)) {
    dir.create(ddir)
    dataUrl <- "https://figshare.com/ndownloader/files/43848153"
    download.file(dataUrl, destfile = dest)
}

Load Seurat object to R environment

bassez <- readRDS(dest)

Explore the dataset and its metadata

# number of genes and cells
dim(bassez)

[1] 20170  9675

# number of cells per patient
table(bassez$sample)


BIOKEY_13_Pre BIOKEY_14_Pre  BIOKEY_5_Pre 
         3369          3185          3121

Metadata content:

bassez@meta.data[1:20, 1:10]

Classify cell types using HiTME

HiTME’s Run.HiTME classification tool is a wrapper of scGate and ProjecTILs to classify cell types in single-cell RNA-seq experiments.

The function takes as input Seurat objects (or list of them). These should be split by sample to avoid batch effects, or split internally in Run.HitME by indicating the parameter split.by.

1st layer

This wrapper firstly runs scGate on TME (Tumor micronenvironment) default models or alternatively on the models provided, resulting in a coarse cell type classification (CD4T, B cell, Dendritic cell…).

2nd layer

Next, HiTME runs ProjecTILs for a finer cell type classification (CD4+ TFH, Tex CD8+, cDC1…) based on the references provided on the cell types classified by scGate that are linked to a respective reference map.

Importantly, HiTME supports analysis using only one of these tools, i.e. only scGate if only coarse cell type classification is intended; or only ProjecTILs, i.e. if a sample is composed exclusively of a high purity coarse cell type subset. This can be achieved by indicating scGate.model = NULL or ref.maps = NULL, respectively.

Additional cell state flavors

Additional signatures not included in the models (e.g. cell cycling, IFN/HSP response…) can also be evaluated using additional.signatures parameter. These are lists of defining genes for these signatures manually costumed or retrieved from SignatuR. As a result, UCell scores for each additional signature will be added as metadata and also summarized when getting a HiT object (see below).

Prepare models and reference maps

Fetch scGate models

require(scGate)

# Fetch all models from scGate
all_models <- get_scGateDB()

# Select models to be used, based on your data specs
HiT_scGate_models <- all_models$human$HiTME

For cell type annotation on whole tumor samples we recommend using the HiTME models from scGate. These contain the following models to classify their respective cell types:

names(HiT_scGate_models)

 [1] "Bcell"         "CD4T"          "CD8T"          "Endothelial"  
 [5] "Epithelial"    "Erythrocyte"   "Fibroblast"    "Mast"         
 [9] "Megakaryocyte" "MoMac"         "Neutrophils"   "NK"           
[13] "panDC"         "PlasmaCell"

These models will be then used for Run.HiTME function (see below). If no model is indicated (scGate.model = "default"), these models will be used by default.

Fetch additional signatures

By default, Run.HiTME will use the Programs group of signatures from SignatuR. Additional custom signatures can be manually added.

require(SignatuR)

species <- "Hs"

additional.signatures <- GetSignature(SignatuR[[species]][["Programs"]])

lapply(additional.signatures, head, n = 3)

$cellCycle.G1S
ENSG00000156802 ENSG00000197299 ENSG00000136492 
        "ATAD2"           "BLM"         "BRIP1" 

$cellCycle.G2M
ENSG00000011426 ENSG00000143401 ENSG00000087586 
         "ANLN"        "ANP32E"         "AURKA" 

$IFN
ENSG00000282851 ENSG00000254838 ENSG00000254781 
        "BISPR"        "GVINP1"        "GVINP2" 

$Tcell.cytotoxicity
[1] "GZMB"  "PRF1"  "FASLG"

$Tcell.exhaustion
[1] "HAVCR2" "CTLA4"  "PDCD1" 

$Tcell.stemness
[1] "TCF7" "SELL" "LEF1"

$M1_macrophage
[1] "CCR7"   "IL2RA"  "IL15RA"

$M2_macrophage
[1] "P2RY13" "LPAR6"  "TGFBR2"

# to add additional signatures

# additional.signatures[['Custom_signature']] <- c('GeneA', 'GeneB', 'GeneC')

Download reference maps

For the finer cell type classification using ProjecTILs we need to download the corresponding reference maps or atlases. These have been previously produced by comprehensive scRNA-seq multi-study integration and curation, and describe reference cell subtypes in a specific TME context. The currently available reference maps for TME context are:

CD4 T cells reference map
CD8 T cells reference map
MoMac (monocytes and macrophages subtypes) reference map
Dendritic cells (DC) reference map

# create directory for reference maps, if not existent
refs_dir <- "refs"
dir.create(refs_dir)

# links to ref maps
ref_links <- c(CD8 = "https://figshare.com/ndownloader/files/41414556", CD4 = "https://figshare.com/ndownloader/files/43794750",
    DC = "https://figshare.com/ndownloader/files/39120752", MoMac = "https://figshare.com/ndownloader/files/42449055")

# download reference maps
lapply(names(ref_links), function(x) {
    download.file(ref_links[[x]], destfile = file.path(refs_dir, paste0(x, ".rds")))
})

Load the reference maps onto the environment

require(ProjecTILs)

ref.maps <- lapply(names(ref_links), function(x) {
    load.reference.map(file.path(refs_dir, paste0(x, ".rds")))
})

[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map Human CD8 TILs"
[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map custom_reference"
[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map DC_human"
[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map MoMac_human"

# give name to the element of the reference list
names(ref.maps) <- names(ref_links)

By default scGate (layer 1) will return the cell ontology ID for each predicted cell type. This ID will be then used to link each coarse cell type with its respective reference map for finer cell type classification using ProjecTILs. Hence, we need to indicate each respective cell ontology ID(s) for each reference map.

If alternative celltype link are used between the coarse and finer cell type classification, this must be specified in Run.HiTME using layer1_link parameter.

# add scGate_link to ref.maps, by default coarse cell type cell ontology ID
layer1.links <- list(CD8 = "CL:0000625", CD4 = "CL:0000624", DC = "CL:0000451", MoMac = "CL:0000576_CL:0000235")
# if multiple coarse cell type are included in a reference map, they should be
# indicated in a vector

for (a in names(ref.maps)) {
    ref.maps[[a]]@misc$layer1_link <- layer1.links[[a]]
}

Before running HiTME, let’s explore the reference maps for the different cell types:

require(ggplot2)
require(patchwork)

refmap_umap <- lapply(names(ref.maps), function(x) {
    DimPlot(ref.maps[[x]], cols = ref.maps[[x]]@misc$atlas.palette, label = T, repel = T) +
        ggtitle(x) + NoLegend()
})

wrap_plots(refmap_umap)

Run HiTME

Now we are set to run Run.HiTME to predict the different coarse and finer cell type found in these breast cancer samples. Ideally, this should be run sample-wise to avoid batch-effect effects on the prediction.

# split Seurat object based on sample
bassez_split <- SplitObject(bassez, split.by = "sample")

We have a list with an element per sample

names(bassez_split)

[1] "BIOKEY_13_Pre" "BIOKEY_14_Pre" "BIOKEY_5_Pre"

Run.HiTME can take a single seurat object or a list of them. Even we can split the object internally using the parameter split.by if the object contains multiple samples. This will return a Seurat object or list of them with new metadata indicating cell type annotation.

bassez_split <- Run.HiTME(bassez_split, scGate.model = HiT_scGate_models, ref.maps = ref.maps)

# Run.HiTME tuning parameters
## do not run
bassez_split <- Run.HiTME(bassez_split,
                          scGate.model = HiT_scGate_models,
                          ref.maps = ref.maps,
                          # already split object
                          split.by = NULL,
                          # if splitting or providing list, whether to return a single merged object
                          remerge = FALSE,
                          # link between scGate and ProjecTILs
                          layer1_link = "CellOntology_ID", 
                          # extra signatures to be computed per celltype
                          additional.signatures = additional.signatures, 
                          # paralelization parameters
                          ncores = 4,
                          progressbar = TRUE
                          )

Get HiT object

Annotated Seurat objects can be summarized into HiT objects using get.HiTObject function. For this function the grouping variable(s) (i.e. cell type classification) group.by resulting from Run.HiTME annotation or additional annotations need to be indicated.

Additional parameters include, useNA parameter which can indicate if undefined/unclassified cells should be considered (NA, by default = FALSE), and clr_zero_impute_perc indicating the % of total cell type counts is added to cell types with 0 when computing CLR (centered-log ratio) transformation (default is 1%).

# create empty list to store HiT objects
hit.list <- list()

# get HiT objects
for(a in names(bassez_split)){
  hit.list[[a]] <- get.HiTObject(bassez_split[[a]],
                                 # multiple grouping variables can be indicated, either as list or vector
                                 group.by = list("layer1" = c("scGate_multi"),
                                       "layer2" = c("functional.cluster"),
                                       "layer1_authors" = c("cellType"),
                                       "layer2_authors" = c("cellSubType")
                                      ),
                                 ## optional parameters
                                 # consider cells not classified (NA)
                                 useNA = T,
                                 # % of total to impute pseudocount for CLR on compositional data
                                 clr_zero_impute_perc = 1
                                 )
}

Alternatively, HiT summarizing object can be obtained directly using Run.HiTME with the parameter return.Seurat = FALSE.

hit.list <- Run.HiTME(bassez_split, ref.maps = ref.maps, remerge = FALSE, return.Seurat = FALSE)

HiT object slots

HiT objects are composed of 4 different slots:

# get one HiT object of the list, in this case corresponds to one sample
hit <- hit.list[[1]]

slotNames(hit)

[1] "metadata"           "predictions"        "composition"       
[4] "aggregated_profile" "version"

Metadata

Seurat object metadata (dataframe): metadata. Including the new metadata added by Run.HiTME and preexisting metadata for each cell.

# see last 10 columns of metadata dataframe
hit@metadata[1:30, (ncol(hit@metadata) - 9):ncol(hit@metadata)]

Predictions

Cell type predictions for each cell in the data set: predictions. This slot is a list where each element is the cell type classification prediction.

names(hit@predictions)

[1] "layer1"         "layer2"         "layer1_authors" "layer2_authors"

# let's inspect the content of layer1, in this case scGate classification
lapply(hit@predictions$layer1, function(x) {
    rownames(x) <- NULL
    head(x, 4)
})

$scGate_multi
[1] "panDC" NA      NA      "panDC"

$additional_signatures
  cellCycle.G1S_UCell cellCycle.G2M_UCell  IFN_UCell Tcell.cytotoxicity_UCell
1          0.00000000          0.03295513 0.05859502                0.0000000
2          0.04766667          0.00000000 0.05725857                0.0000000
3          0.00000000          0.03359615 0.05236449                0.2883333
4          0.00000000          0.03200000 0.12772897                0.0000000
  Tcell.exhaustion_UCell Tcell.stemness_UCell M1_macrophage_UCell
1                      0                    0          0.04839506
2                      0                    0          0.05733333
3                      0                    0          0.03543210
4                      0                    0          0.15191975
  M2_macrophage_UCell
1          0.05982171
2          0.04720930
3          0.00000000
4          0.11909302

$scGate_is.pure
  is.pure_Bcell is.pure_CD4T is.pure_CD8T is.pure_Endothelial
1        Impure       Impure       Impure              Impure
2        Impure       Impure       Impure              Impure
3        Impure       Impure       Impure              Impure
4        Impure       Impure       Impure              Impure
  is.pure_Epithelial is.pure_Erythrocyte is.pure_Fibroblast is.pure_Mast
1             Impure              Impure             Impure       Impure
2             Impure              Impure             Impure       Impure
3             Impure              Impure             Impure       Impure
4             Impure              Impure             Impure       Impure
  is.pure_Megakaryocyte is.pure_MoMac is.pure_Neutrophils is.pure_NK
1                Impure        Impure              Impure     Impure
2                Impure        Impure              Impure     Impure
3                Impure        Impure              Impure     Impure
4                Impure        Impure              Impure     Impure
  is.pure_panDC is.pure_PlasmaCell
1          Pure             Impure
2        Impure             Impure
3        Impure             Impure
4          Pure             Impure

$UCell_scores
  Immune_UCell Lymphoid_UCell PanBcell_UCell Bcell_UCell APC_UCell
1    0.3592000              0              0           0  0.968625
2    0.5113333              0              0           0  0.000000
  Epithelial_UCell Stromal_UCell Myeloid_UCell MoMacDC_UCell Neutrophils_UCell
1                0     0.1205238             0     0.1448571                 0
2                0     0.0000000             0     0.0000000                 0
  Tcell_UCell NK_UCell Endothelial_UCell Talphabeta_UCell Plasma_cell_UCell
1       0.000   0.0000                 0        0.0000000                 0
2       0.409   0.0539                 0        0.2386667                 0
  CD4Tcell_UCell Treg_UCell CD8T_UCell Tgammadelta_UCell TAM_UCell Tgd1_UCell
1      0.0000000          0          0         0.0000000         0  0.0000000
2      0.2736667          0          0         0.3404667         0  0.8481667
  Tgd2_UCell CD4T_UCell Fibroblast_UCell Alveolar_type1_UCell
1  0.0000000          0        0.1205238                    0
2  0.3573333          0        0.0000000                    0
  Alveolar_type2_UCell Erythrocyte_UCell Mast_UCell Megakaryocyte_UCell
1                    0                 0          0                   0
2                    0                 0          0                   0
  Macrophage_UCell Monocyte_UCell pDC_UCell cDC1_UCell cDC2_UCell cDC2b_UCell
1                0              0 0.2084167          0  0.6062222   0.2084167
2                0              0 0.0000000          0  0.0000000   0.0000000
  DC3_UCell Monocyte_not_cDC2_UCell
1 0.1205238                       0
2 0.0000000                       0
 [ reached 'max' / getOption("max.print") -- omitted 2 rows ]

# and also for coarse classification by authors (cellType)
lapply(hit@predictions$layer1_authors, function(x) {
    rownames(x) <- NULL
    head(x, 4)
})

$cellType
[1] "Myeloid_cell" "T_cell"       "pDC"          "Myeloid_cell"

Confusion matrix

We can visualize the agreement between HiTME celltype prediction and authors manual annotation in Bassez et al. (2021) using a confusion matrix with HiTME function plot.confusion.matrix.

# use all 3 samples (whole list of HiT objects)
plot.confusion.matrix(hit.list,
                      var.1 = "cellType",
                      var.2 = "scGate_multi",
                      # use relative proportion of cells
                      relative = TRUE,
                      # produce confusion matrix plot
                      type = "tile")

We cal also indicate type = "barplot" to visualize the agreement using a bar plot:

# use relative proportion of cells
plot.confusion.matrix(hit.list, var.1 = "cellType", var.2 = "scGate_multi", relative = TRUE,
    type = "barplot")

Composition

Cell type composition for each layer of cell type prediction: composition. A list including:

Cell type counts (cell_counts)
Cell type relative proportions (freq)
Centered log ratio (CLR)-transformed cell type counts (freq_clr)

# let's inspect the compositional data for HiTME layer 1 (scGate)
hit@composition$layer1

$cell_counts
  Bcell CD4T CD8T Endothelial Epithelial Fibroblast Mast MoMac NK panDC
1   115  819  669         242         19        841    6    47 58    87
  PlasmaCell  NA
1         47 419

$freq
     Bcell     CD4T     CD8T Endothelial Epithelial Fibroblast      Mast
1 3.413476 24.30988 19.85752     7.18314  0.5639656    24.9629 0.1780944
     MoMac       NK    panDC PlasmaCell       NA
1 1.395073 1.721579 2.582369   1.395073 12.43692

$freq_clr
       Bcell     CD4T     CD8T Endothelial Epithelial Fibroblast      Mast
1 -0.2107966 1.535736 1.342256   0.4066168  -1.248235   1.561209 -1.531561
       MoMac         NK      panDC PlasmaCell        NA
1 -0.8220456 -0.6942469 -0.4194349 -0.8220456 0.9025474

Cell type frequency plots

We can visualize the composition data of all 3 samples using the HiTME function plot.celltype.freq. The parameter by.x defines which variables, sample or cell type, are included in the x-axis the plot types, rendering a barplot or boxplot, respectively.

# use all 3 samples (whole list of HiT objects)

# by sample, barplot
plot.celltype.freq(hit.list, group.by = "scGate_multi", by.x = "sample")

# by sample, boxplot
plot.celltype.freq(hit.list, group.by = "scGate_multi", by.x = "celltype")

Aggregated profile

The aggregated profile slot of the HiT object is composed of a list of 2 elements:

Pseudobulk matrix: Aggregated counts per cell type corresponding of each classification.
Averaged signature scores: Average of signature score for indicated additional signatures for each cell type classified.

# content of aggregated_profile slot of HiT object
names(hit@aggregated_profile)

[1] "Pseudobulk" "Signatures"

# first entries of each pseudobulk matrix gene expression for each
# classification (group.by)
lapply(hit@aggregated_profile$Pseudobulk, function(x) {
    x[1:3, 1:3]
})

$layer1
3 x 3 sparse Matrix of class "dgCMatrix"
         Bcell CD4T CD8T
A1BG         7   57   56
A1BG-AS1     5   13    9
A2M          .    8    6

$layer2
3 x 3 sparse Matrix of class "dgCMatrix"
         AS-DC CD4.CTL-EOMES CD4.CTL-GNLY
A1BG         1             5            7
A1BG-AS1     .             1            .
A2M          .             .            .

$layer1_authors
3 x 3 sparse Matrix of class "dgCMatrix"
         B-cell Cancer-cell Endothelial-cell
A1BG         34           3               33
A1BG-AS1      6           .                1
A2M           4          18             1905

$layer2_authors
3 x 3 sparse Matrix of class "dgCMatrix"
         ASDC C1-CD14 C10-LYVE1
A1BG        1       .         .
A1BG-AS1    .       .         .
A2M         .       .        19

# first entries of each average additional signature UCell score for each
# classification (group.by)
lapply(hit@aggregated_profile$Signatures, function(x) {
    x[1:3, 1:4]
})

$layer1
# A tibble: 3 × 4
  scGate_multi cellCycle.G1S_UCell cellCycle.G2M_UCell IFN_UCell
  <chr>                      <dbl>               <dbl>     <dbl>
1 Bcell                     0.0282              0.0286    0.0657
2 CD4T                      0.0248              0.0287    0.0670
3 CD8T                      0.0234              0.0290    0.0688

$layer2
# A tibble: 3 × 4
  functional.cluster cellCycle.G1S_UCell cellCycle.G2M_UCell IFN_UCell
  <chr>                            <dbl>               <dbl>     <dbl>
1 AS-DC                           0.0278              0.0309    0.0937
2 CD4.CTL_EOMES                   0.0251              0.0336    0.0640
3 CD4.CTL_GNLY                    0.0222              0.0268    0.0700

$layer1_authors
# A tibble: 3 × 4
  cellType         cellCycle.G1S_UCell cellCycle.G2M_UCell IFN_UCell
  <chr>                          <dbl>               <dbl>     <dbl>
1 B_cell                        0.0246              0.0258    0.0696
2 Cancer_cell                   0.0323              0.0360    0.0714
3 Endothelial_cell              0.0228              0.0264    0.0870

$layer2_authors
# A tibble: 3 × 4
  cellSubType cellCycle.G1S_UCell cellCycle.G2M_UCell IFN_UCell
  <chr>                     <dbl>               <dbl>     <dbl>
1 ASDC                     0.0278              0.0309    0.0937
2 C10_LYVE1                0.0234              0.0266    0.0828
3 C1_CD14                  0                   0.0328    0.0919

Overall cell type classification evaluation

We can observe a high degree of agreement between HiTME coarse cell type prediction and manual annotation by authors (Bassez et al., 2021) on these 3 samples:

B cells: 70% of agreement (considering that authors do not classify plasma cells).
Endothelial cells: 90% of agreement.
Fibroblasts: 90% of agreement.
Mast cells: 100% of agreement.
Myeloid cells (MoMac): 40% of agreement. Here HiTME shows apparently discrepancies with authors classification, and classify 30% of author’s Myeloid cells as panDC. However, most DC subtypes are indeed myeloid, and authors only provide classification for plasmacytoid DC.
Dendritic cells: 90% of agreement.
T cells: 80% of agreement, considering that authors do not split between CD4+ and CD8+ T cells.

Finally, HiTME still do not support cancer cell classification. Most of the cancer cells regarded by authors (confirmed by copy number-unstability) are undefined by HiTME as expected (50%). However, 40% of them are classified as epithelial cells, as these cancerous cells come from a breast cancer.

Regarding the cell type composition of these tumor samples, we can see that all 3 samples are mostly composed of T cells, with larger proportion of CD4+ T cells vs CD8+ counterparts. The second most common cell type are fibroblasts, albeit with high variability, followed by B cells. The number of Monocytes, macrophages and dendritic cells is low in these samples.

Visualization of cell type classification

We can also integrate the 3 samples using STACAS to minimize possible batch effects and evaluate cell type clustering on a low-dimensionality space.

# integrate with STACAS
require(STACAS)

# set PCA dimensions to use
dims <- 30

bassez_integrated <- bassez_split %>%
    Run.STACAS(dims = dims) %>%
    # Run UMAP on integrated PCA space
RunUMAP(dims = 1:dims)

Show UMAP comparing HiTME and author’s coarse classification :

require(ggplot2)

# show HiTME and authors coarse cell type classification
grouping <- c("scGate_multi", "cellType")
umap.list <- lapply(grouping, function(g) {
    DimPlot(bassez_integrated, group.by = g, label = T, repel = T) + guides(color = guide_legend(ncol = 3,
        override.aes = list(size = 3))) + theme(legend.position = "bottom")
})

wrap_plots(umap.list)

With this visualization we can observe firstly that authors classified few B cells that are in fact clustering with other cell types (Fibroblasts, endothelial cells, and myeloid cells). Moreover, the cell type with highest discrepancy between HiTME and authors, MoMac vs myeloid cells, seem to be actually forming two different subclusters. Actually, authors only classify a distinct cluster of plasmacytoid dendritic cells (pDC), but do not show classification for myeloid dentritic cells.

We can explore this further by checking the markers on these myeloid cells.

DefaultAssay(bassez_integrated) <- "RNA"

# keep only MoMac and panDC classification to compare them
bassez_momac <- bassez_integrated[, bassez_integrated$cellType == "Myeloid_cell" &
    bassez_integrated$scGate_multi %in% c("MoMac", "panDC")]

momacDC_markers <- c("S100A9", "S100A8", "FCN1", "VCAN", "C1QA", "C1QB", "APOE",
    "APOC1", "C1QC", "CD68", "FCER1G", "CLEC9A", "XCR1", "IRF8", "CLEC10A", "FCER1A",
    "CD1C")

# keep only genes expressed
sc_names = rownames(bassez_momac)[rowSums(bassez_momac) > 0]

momacDC_markers <- intersect(momacDC_markers, sc_names)

VlnPlot(bassez_momac, features = momacDC_markers, group.by = "scGate_multi", stack = TRUE,
    flip = TRUE, ncol = 2) + NoLegend() + ggtitle("Author's Myeloid cells") + xlab("HiTME classification")

With this plot we can see that “Myeloid cells” according to authors, not only seem to form two independent clusters on the UMAP but also have an striking different expression of MoMac and dendritic cells markers in each cell type.

Visualization of finer classification

As we have seen, HiTME and authors classification show high agreement on overall T cell classification (~80%). However, authors do not tell apart CD4+ and CD8+ T cells, and do not take into account NK cells.

Let’s inspect how HiTME and author’s classify T cells at finer level

require(dplyr)

# filter Hit metadata for only T cells and remove finer classification NA
hit_Tcell <- lapply(hit.list, function(x) {
    x@metadata <- x@metadata %>%
        filter(scGate_multi %in% c("CD4T", "CD8T") & !is.na(functional.cluster) &
            !is.na(cellSubType))
    return(x)
})

Confusion matrix

plot.confusion.matrix(hit_Tcell,
                      var.1 = "cellSubType",
                      var.2 = "functional.cluster",
                      # use relative proportion of cells
                      relative = TRUE,
                      # produce confusion matrix plot
                      type = "tile")

Integrated UMAP

tcell_integrated <- bassez_split %>%
    # filter for author T cells
lapply(., function(x) {
    x[, x$cellType == "T_cell"]
}) %>%
    Run.STACAS(dims = dims) %>%
    # Run UMAP on integrated PCA space
RunUMAP(dims = 1:dims)

# show both finer cell type classification and also coarse of HiTME
grouping <- c("scGate_multi", "functional.cluster", "cellSubType")
umap.list <- lapply(grouping, function(g) {
    DimPlot(tcell_integrated, group.by = g, label = T, repel = T) + guides(color = guide_legend(ncol = 3,
        override.aes = list(size = 3))) + theme(legend.position = "bottom")
})

wrap_plots(umap.list, ncol = 2, nrow = 2)

Overall, for finer cell types is harder to evaluate the classification performance, given the lack of consensus on the subtype names and definition. However, HiTME and authors seem to largely agree on most subtypes.

References

Andreatta, Massimo, Ariel J. Berenstein, and Santiago J. Carmona. 2022. “scGate: Marker-Based Purification of Cell Types from Heterogeneous Single-Cell RNA-Seq Datasets.” Bioinformatics 38 (April): 2642–44. https://doi.org/10.1093/BIOINFORMATICS/BTAC141.

Andreatta, Massimo, Jesus Corria-Osorio, Sören Müller, Rafael Cubas, George Coukos, and Santiago J. Carmona. 2021. “Interpretation of t Cell States from Single-Cell Transcriptomics Data Using Reference Atlases.” Nature Communications 2021 12:1 12 (May): 1–19. https://doi.org/10.1038/s41467-021-23324-4.

Andreatta, Massimo, Léonard Hérault, Paul Gueguen, David Gfeller, Ariel J Berenstein, and Santiago J Carmona. 2023. “Semi-Supervised Integration of Single-Cell Transcriptomics Data.” bioRxiv, July, 2023.07.07.548105. https://doi.org/10.1101/2023.07.07.548105.

Bassez, Ayse, Hanne Vos, Laurien Van Dyck, Giuseppe Floris, Ingrid Arijs, Christine Desmedt, Bram Boeckx, et al. 2021. “A Single-Cell Map of Intratumoral Changes During Anti-PD1 Treatment of Patients with Breast Cancer.” Nature Medicine 2021 27:5 27 (May): 820–32. https://doi.org/10.1038/s41591-021-01323-8.

Cell type classification in tumor microenvironment using HiTME